Distributed-memory Hierarchical Interpolative Factorization

The hierarchical interpolative factorization (HIF) offers an efficient way for solving or preconditioning elliptic partial differential equations. By exploiting locality and low-rank properties of the operators, the HIF achieves quasi-linear complexity for factorizing the discrete positive definite elliptic operator and linear complexity for solving the associated linear system. In this paper, the distributed-memory HIF (DHIF) is introduced as a parallel and distributed-memory implementation of the HIF. The DHIF organizes the processes in a hierarchical structure and keep the communication as local as possible. The computation complexity is $O\left(\frac{N\log N}{P}\right)$ and $O\left(\frac{N}{P}\right)$ for constructing and applying the DHIF, respectively, where $N$ is the size of the problem and $P$ is the number of processes. The communication complexity is $O\left(\sqrt{P}\log^3 P\right)\alpha + O\left(\frac{N^{2/3}}{\sqrt{P}}\right)\beta$ where $\alpha$ is the latency and $\beta$ is the inverse bandwidth. Extensive numerical examples are performed on the NERSC Edison system with up to 8192 processes. The numerical results agree with the complexity analysis and demonstrate the efficiency and scalability of the DHIF

Bold Diagrammatic Monte Carlo in the Lens of Stochastic Iterative Methods

This work aims at understanding of bold diagrammatic Monte Carlo (BDMC) methods for stochastic summation of Feynman diagrams from the angle of stochastic iterative methods. The convergence enhancement trick of the BDMC is investigated from the analysis of condition number and convergence of the stochastic iterative methods. Numerical experiments are carried out for model systems to compare the BDMC with related stochastic iterative approaches

Interpolative Butterfly Factorization

This paper introduces the interpolative butterfly factorization for nearly optimal implementation of several transforms in harmonic analysis, when their explicit formulas satisfy certain analytic properties and the matrix representations of these transforms satisfy a complementary low-rank property. A preliminary interpolative butterfly factorization is constructed based on interpolative low-rank approximations of the complementary low-rank matrix. A novel sweeping matrix compression technique further compresses the preliminary interpolative butterfly factorization via a sequence of structure-preserving low-rank approximations. The sweeping procedure propagates the low-rank property among neighboring matrix factors to compress dense submatrices in the preliminary butterfly factorization to obtain an optimal one in the butterfly scheme. For an $N\times N$ matrix, it takes $O(N\log N)$ operations and complexity to construct the factorization as a product of $O(\log N)$ sparse matrices, each with $O(N)$ nonzero entries. Hence, it can be applied rapidly in $O(N\log N)$ operations. Numerical results are provided to demonstrate the effectiveness of this algorithm

Globally Constructed Adaptive Local Basis Set for Spectral Projectors of Second Order Differential Operators

Spectral projectors of second order differential operators play an important role in quantum physics and other scientific and engineering applications. In order to resolve local features and to obtain converged results, typically the number of degrees of freedom needed is much larger than the rank of the spectral projector. This leads to significant cost in terms of both computation and storage. In this paper, we develop a method to construct a basis set that is adaptive to the given differential operator. The basis set is systematically improvable, and the local features of the projector is built into the basis set. As a result the required number of degrees of freedom is only a small constant times the rank of the projector. The construction of the basis set uses a randomized procedure, and only requires applying the differential operator to a small number of vectors on the global domain, while each basis function itself is supported on strictly local domains and is discontinuous across the global domain. The spectral projector on the global domain is systematically approximated from such a basis set using the discontinuous Galerkin (DG) method. The global construction procedure is very flexible, and allows a local basis set to be consistently constructed even if the operator contains a nonlocal potential term. We verify the effectiveness of the globally constructed adaptive local basis set using one-, two- and three-dimensional linear problems with local potentials, as well as a one dimensional nonlinear problem with nonlocal potentials resembling the Hartree-Fock problem in quantum physics

Triangularized Orthogonalization-free Method for Solving Extreme Eigenvalue Problems

A novel orthogonalization-free method together with two specific algorithms are proposed to solve extreme eigenvalue problems. On top of gradient-based algorithms, the proposed algorithms modify the multi-column gradient such that earlier columns are decoupled from later ones. Global convergence to eigenvectors instead of eigenspace is guaranteed almost surely. Locally, algorithms converge linearly with convergence rate depending on eigengaps. Momentum acceleration, exact linesearch, and column locking are incorporated to further accelerate both algorithms and reduce their computational costs. We demonstrate the efficiency of both algorithms on several random matrices with different spectrum distribution and matrices from computational chemistry

Variational training of neural network approximations of solution maps for physical models

A novel solve-training framework is proposed to train neural network in representing low dimensional solution maps of physical models. Solve-training framework uses the neural network as the ansatz of the solution map and train the network variationally via loss functions from the underlying physical models. Solve-training framework avoids expensive data preparation in the traditional supervised training procedure, which prepares labels for input data, and still achieves effective representation of the solution map adapted to the input data distribution. The efficiency of solve-training framework is demonstrated through obtaining solutions maps for linear and nonlinear elliptic equations, and maps from potentials to ground states of linear and nonlinear Schr\"odinger equations

Coordinate-wise descent methods for leading eigenvalue problem

Leading eigenvalue problems for large scale matrices arise in many applications. Coordinate-wise descent methods are considered in this work for such problems based on a reformulation of the leading eigenvalue problem as a non-convex optimization problem. The convergence of several coordinate-wise methods is analyzed and compared. Numerical examples of applications to quantum many-body problems demonstrate the efficiency and provide benchmarks of the proposed coordinate-wise descent methods

A Multiscale Butterfly Algorithm for Multidimensional Fourier Integral Operators

This paper presents an efficient multiscale butterfly algorithm for computing Fourier integral operators (FIOs) of the form (\mathcal{L} f)(x) = \int_{R^d}a(x,\xi) e^{2\pi \i \Phi(x,\xi)}\hat{f}(\xi) d\xi, where $\Phi(x,\xi)$ is a phase function, $a(x,\xi)$ is an amplitude function, and $f(x)$ is a given input. The frequency domain is hierarchically decomposed into a union of Cartesian coronas. The integral kernel a(x,\xi) e^{2\pi \i \Phi(x,\xi)} in each corona satisfies a special low-rank property that enables the application of a butterfly algorithm on the Cartesian phase-space grid. This leads to an algorithm with quasi-linear operation complexity and linear memory complexity. Different from previous butterfly methods for the FIOs, this new approach is simple and reduces the computational cost by avoiding extra coordinate transformations. Numerical examples in two and three dimensions are provided to demonstrate the practical advantages of the new algorithm

Multidimensional Butterfly Factorization

This paper introduces the multidimensional butterfly factorization as a data-sparse representation of multidimensional kernel matrices that satisfy the complementary low-rank property. This factorization approximates such a kernel matrix of size $N\times N$ with a product of \O(\log N) sparse matrices, each of which contains \O(N) nonzero entries. We also propose efficient algorithms for constructing this factorization when either (i) a fast algorithm for applying the kernel matrix and its adjoint is available or (ii) every entry of the kernel matrix can be evaluated in \O(1) operations. For the kernel matrices of multidimensional Fourier integral operators, for which the complementary low-rank property is not satisfied due to a singularity at the origin, we extend this factorization by combining it with either a polar coordinate transformation or a multiscale decomposition of the integration domain to overcome the singularity. Numerical results are provided to demonstrate the efficiency of the proposed algorithms

On the numerical rank of radial basis function kernels in high dimension

Low-rank approximations are popular methods to reduce the high computational cost of algorithms involving large-scale kernel matrices. The success of low-rank methods hinges on the matrix rank of the kernel matrix, and in practice, these methods are effective even for high-dimensional datasets. Their practical success motivates our analysis of the function rank, an upper bound of the matrix rank. In this paper, we consider radial basis functions (RBF), approximate the RBF kernel with a low-rank representation that is a finite sum of separate products and provide explicit upper bounds on the function rank and the $L_\infty$ error for such approximations. Our three main results are as follows. First, for a fixed precision, the function rank of RBFs, in the worst case, grows polynomially with the data dimension. Second, precise error bounds for the low-rank approximations in the $L_\infty$ norm are derived in terms of the function smoothness and the domain diameters. Finally, a group pattern in the magnitude of singular values for RBF kernel matrices is observed and analyzed, and is explained by a grouping of the expansion terms in the kernel's low-rank representation. Empirical results verify the theoretical results